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Abstract 

We propose a mathematical model which suggests that the two main 
geological observations about shingle beaches, i.e. the emergence of pre- 
dominant pebble size ratios and strong segregation by size are interrelated. 
Our model is a based on a system of ODEs called the 60a; equations, 
describing the evolution of pebble ratios. We derive these ODEs as a 
heuristic approximation of Bloore's PDE describing coUisional abrasion 
and verify them by simple experiments and by direct simulation of the 
PDE. While representing a radical simplification of the latter, our system 
admits the inclusion of additional terms related to frictional abrasion. We 
show that nontrivial attractors (corresponding to predominant pebble size 
ratios) only exist in the presence of friction. By interpreting our equa- 
tions as a Markov process, we illustrate by direct simulation that these 
attractors may only be stabilized by the ongoing segregation process. 

1 Introduction 



The shape of pebbles has been a matter of discussion since at least the time 
of Aristotle [TH]- In general, the central question is whether particular pebble 
shapes emerge from the abrasion and transport processes. Aristotle himself 
claimed that spherical shapes dominate. However, as Aristotle also observed, 
abrasion is a complex interaction between the abraded particle and the abrading 
environment represented by 'other objects' (i.e. other pebbles) where not only 
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local properties (e.g. curvatures) but also semi-global effects due to particle 
shape as well as global effects due to particle transport play important role. 
In this process pebbles mutually abrade each other, defining the time evolu- 
tion both for the abraded pebble and the abrading environment represented by 
other particles subject to particle transport. In the simplest approach, one ne- 
glects the latter effects and regards the abrasion of a single pebble in a constant 
environment. Aristotle's model for individual abrasion may be translated as 

« - fiR) , (1) 

where v is the speed with which the pebble's surface moves along the inward 
normal, R is the radial distance from the center of gravity of the pebble being 
abraded, / is a monotonically increasing function of R only, and in particular, 
/ is independent of time. Note that ([T]) is a partial integro-differential equation, 
since the location of the center of gravity is determined by time-dependent 
integrals. The modern , PDE-based, theory of individual abrasion appears 
to start with the pioneering work, both experimental and theoretical, of Lord 
Rayleigh [24l |25l |26] . For some earlier work see [30l [8] . Rayleigh mainly 
considered axisymmetric pebbles and he observed that the ultimate shape was 
not necessarily, indeed seldom, spherical. He found that some pebble's shapes 
are far from ellipsoidal, being much more discoid in shape. He asserted that 
abrasion cannot be a simple function of the Gauss (or specific) curvature K . 

On the mathematical side, Firey [10] initiated a rigorous study by adopting 
a PDF model rejected by Rayleigh in which the shape evolved according to what 
Andrews |T] calls the How by the Gauss curvature that is, he studied the PDF 

v = cK, (2) 

where c is a constant and proved that all convex shapes ultimately converge 
to the sphere under the action of ([2]). Note that the word flow is being used 
in the sense in dynamical systems theory and it should not be confused with 
physical fluid flow. Later, this proof was substantially amplified by Andrews 
([T]). Recently, Durian [5] investigated the statistical distribution of Gaussian 
curvature on pebble shapes. 

The physical assumption underlying Firey's model is that the abraded par- 
ticle (pebble) undergoes a series of small collisions with a very large, smooth 
abrader, and this might be the case when pebbles are carried by a fast river 
and colhde repeatedly with the riverbed, a process called bedload transport. 
This concept of collisional flows was substantially generalized by Bloore who 
studied abrasion by spherical abraders of radius r (Firey's model corresponds 
to r — > oo). With brilliant intuition Bloore [3] arrived at the PDF 

v = l + 2bH + cK , (3) 

where b and c represent, as Bloore describes, the average, or 'effective' values 
for r and r^, respectively and H is the Mean Curvature. Bloore also showed 

^Son of the Nobel Prize winner Lord Rayleigh 



2 



that a spherical shape of radius R abraded by an abrader of radius r is stable 
under the action of ([3]) if 

-<3. (4) 

r 

Apparently, Q consists of three terms: the constant first term which corre- 
sponds to the so-called Eikonal equation, the Mean Curvature term 2bH and 
the Gaussian term cK. The latter was studied by Firey, while for the Mean 
Curvature Flow Huisken fT4| showed that it also converges to the sphere. In the 
case of the Eikonal we know that the outward flow (which can be imagined also 
as time-reversed abrasion) converges to the sphere, so, simply put, the question 
is whether the three terms can 'balance' i.e. homothetic solutions can exist. 

The rigorous foundation of the collisional model requires a probabilistic ap- 
proach as it was already conjectured by Hilbert [13]. Based on results from 
Stochastic and Integral Geometry 22. one can not only justify Bloore's equa- 
tions but also identify the physical meaning of the coefficients 6, c in ^ [S^ : 

where M is the integrated mean curvature and A the area of the abrading parti- 
cles, assumed convex and identical. In case of spherical abraders this is identical 
with Bloore's interpretation of the coefhcients. By Minkowski's inequality for 
convex bodies |27j the quantities b and c may not be given arbitrarily, rather, 
based on ([5]) they must satisfy 

b^>c. (6) 

The main question is now whether collisional flows defined in ([3]) , ([5]) can explain 
the observation of spherical and non-spherical pebbles shapes, i.e. whether ([3]), 
([5]) admits homothetic solutions and what they look like. 

Independently of the mathematical work which we have briefly reviewed, 
there has been considerable geological interest in pebble shape which partly 
focused on individual pebble shapes, partly on global understanding of pebble 
transport and spatial distribution on large shingle beaches. An extensive list of 
references about individual pebble shape up to 1970 is contained in [8]. 

One of the most remarkable accounts of pebble shape and distribution is 
provided by Carr [6] based on the measurement of approximately a hundred 
thousand pebbles on Chesil Beach, Dorset, England. In summarizing his re- 
sults, Carr first notes that a pebble appears to retain its proportions irrespective 
of size and angularity. Apparently, the same phenomenon was described by Lan- 
don [20j claiming that coastal pebbles evolve towards what he calls equilibrium 
shapes.The subject was taken up in the 1990's in Nature O [22l following a 
note by Wald [S^ claiming that pebbles with three planes of symmetry typically 
had the form of flattened ovoids in proportion 7:6:3. In the absence of pebble 
transport, these observations would imply the existence of stable (attracting) 
homothetic (self-similar) solutions of the governing equations. As we will point 
out, in the absence of transport such solutions do not exist. 
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However, Carr also provides mean values and sample variations for maximal 
pebble size along lines orthogonal to the beach. These plots reveal pronounced 
segregation by maximal size, i.e. on shingle beaches pebbles of roughly similar 
maximal sizes appear to be spatially close to each other. This phenomenon, also 
confirmed by [3], [II],[I1], [IH],[in], HH] is closely related to the global transport 
of pebbles by waves [H], [7]. Grading not only by size but also by shape has 
been also reported, the most notable description is due to Bluck 0, reporting on 
what he calls equilibrium states of pebble collections corresponding to stationary 
shape distributions (cf also ^U\, [53], [35]). 

Apparently, the distribution of pebble shapes and size are controlled by 
two fundamental geological processes: particle abrasion and particle transport. 
Which of these processes dominates may depend on the geological location and 
also on time, however, geologists appear to agree that the role of transport may 
not be neglected. Another general observation is that this process approaches 
equilibria in the form of stationary pebble shape distributions with dominant 
peaks. A detailed account of the interaction of abrasion and transport is given 
by Landon [50] investigating the beaches on the west shore of Lake Michigan. 
He attributes shape variation to a mixture of abrasion and transport. Kuenen 
[T^ discusses Landon's observations, however disagrees with the conclusions and 
attributes shape variation primarily to transport. Carr [5] observes dominant 
shape ratios emerging as a result of abrasion and size grading while Bluck [S] 
describes beaches in South Wales where equilibrium distributions of shape are 
reached primarily by transport and abrasion plays a minor role. 

These observations call for a mathematical model where both abrasion and 
transport are included and one may study their complex interaction in the 
entire range from pure abrasion to pure segregation. Our goal in this paper is 
to make one step towards such a general model which is based on the classical 
theory of particle abrasion but also admits the inclusion of global transport. 
We derive a system of representative ODE's called the box equations which 
provide a huge simplification of Bloore's PDE's but which nevertheless allow us 
to study the spatio-temporal process. In this framework we can show that in 
accord with the observations of Landon, Carr and Bluck, in the absence of global 
transport and segregation, stationary distributions with sharp peaks (centered 
around homothetic solutions of the deterministic equations) may not emerge in a 
stable manner. After verifying our analytical results both by direct simulation 
of Bloore's PDE and by simple experiments we also show that abrasion and 
size grading alone (most clearly stated by Carr) might be sufficient to produce 
persistent dominant basic proportions. Strong shape grading (as described by 
Bluck) enhances this process even further. Our model creates a framework 
where both extreme cases (pure abrasion in the absence of transport and pure 
transport in the absence of abrasion), as well as their combination may be 
studied. 
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2 The CoUisional Box Equations 



Motivated partly by the above described geological observations, instead of the 
numerical study of (jSj we propose representing shapes by orthogonal bounding 
boxes (defined by 6 orthogonal bounding planes) with sizes 2ui < 2u2 < 2w3, so 
that the attrition speed v is replaced by dui/dt — Ui, the local quantities K,H 
we obtain from the ellipsoid which fits inside the box: 

K,^^, H,^U^+'^], ^^J^k. (7) 
ujul 2 \^u2 ulj 

Since we are interested in the existence of self-similar solutions, we introduce 
the box ratios yi — ui/u^, j/2 — ^2/^3, we denote the origin (j/i, 1/2) = (0, 0) by 
O the sphere at (yi, y2) = (l,l)byS' and to simplify notation we write 2/3 = U3. 
This approach has the advantage that most of the available geological data is 
expressed in terms of these variables, in fact, the representation in the (2/1,2/2) 
plane is often referred to as the Zingg Triangle [37] . Now the PDE ([3]) is reduced 
to the system of O.D.E,'s in the (2/1,2/2,2/3) system: 



pE pM pG 

m = i^^(2/i,2/2,y3,&,c) = ^ + 26^+c^ (8) 

2/3 2/3 2/3 

-m = F3(2/i,2/2,2/3,&,c) = IH % 2^ + ' (9) 

2/3 2/12/2 2/5 2/12/2 



where 

Ff^y,-1, Ft'=^—^ = *,J = l,2;z^j. (10) 

22/^ 2/^2/^ 



As can be seen from ([H])-®, the full flow is three dimensional, however, the 
component flows (Eikonal, Mean curvature, Gaussian) have autonomous planar 
components (illustrated in figure [Ij: 

,E_Fl_ y2^ ,M _Fi^_ 2/1(2/1-1) ,g_EL- vM-^) 
2/1-1' Ff 2/2(2/?-l)' F^ vM-l)' 

(11) 

where ()' = d/dyi. By introducing the vector notation y — [2/1, 2/2, 2/3]^, F = 
[i^i,F2,^3]'^, dS])-® can be rewritten as 

y = F(y,6,c). (12) 

Our main concern is the detection of fixed points and self-similar (homothetic) 
solutions. According to © 2/3 < 0, so ([S])-® does not have fixed points. Ho- 
mothetic solutions satisfying 

2/1 (0 = 2/2(0 = (13) 
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may exist, in fact spheres at S are a homothetic solution of (IH])-®. Apart 
from spheres, ([8])- ([9]) will only admit homothetic solutions if (H)) is autonomous 
and has fixed points; we explore this case in subsection 12.21 If ([S]) is not au- 
tonomous, then homothetic solutions do not exist, nevertheless one still may 
look for solutions of 

yi{to) = hito) = (14) 

as 'temporary' homothetic solutions, we explore this in subsection l2.1l Whether 
permanent or temporary, the existence of homothetic solutions requires that the 
(yij2/2) projection of the flow should be balanced and this observation admits 
the formulation of a necessary condition. By considering that the flow is a linear 
combination of three components (jlll) two of which converge to sphere, we can 
observe that if either or (|14p are to be satisfied then the slope of the Eikonal 
should be 'bracketed' by the two other slopes: 

-yr)<o. (15) 

Based on (jlip we can readily verify that in the interior j/i < ?/2 < 1 of the Zingg 
triangle this is indeed the case since 

vr < < 2^2^ < I ' (16) 

i.e. the necessary condition (fT5|) is always fulfilled, cf. Figure |4)/a. 




Figure 1: 2D component flows: Eikonal, Mean Curvature and Gaussian flows 

Box flows represent only a heuristic approximation of the PDF and we may 
use several alternative definitions of the bounding box when measuring the 
box dimensions. Selection of special approximations based on extremal proper- 
ties has been studied and a famous example is the unique John ellipsoid with 
maximal volume [T5]. Box equations are the basic tool enabling us to follow 
Aristotle's ideas and to include semi-global and global effects in our model. On 
the other hand, box equations are Aristotelean also in the local sense since in 
the box approximations of the Bloore flows ([3]) attrition speeds dui/dt — Ui 
are a monotonically growing function of the distance ui. When deriving the 
box ode's from the PDF the abrading environment (represented by the 
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coefficients b, c) can also be approximated by its box dimensions vi < V2 < 
and box ratios zi — vi/v^, Z2 — vijv-i and maximal box size Z3 = W3. According 
to ([5]) the coefficients 6, c depend on the surface integrals M, A, and the latter 
we derive directly from the abrading orthogonal box, with ratios Zi. 

M = 27rz3(zi + Z2 + 1), ^ = 8zf(zi2;2 + 21 + ^2). (17) 

The same quantities can be expressed for the unit cube as M\ — Gn, Ai — 24, 
so we have 

U \ Z1+Z2 + I f, , . A 2 ^1 + Z2 + Z1Z2 2 rc 

b{.) = — ^z, = .3/., c(z) = - = Z3 = .3/.- 

. 

We remark that ([6]) remains valid also in the box approximations. While the 
PDE model ([3]) is describing the collisional abrasion of an individual pebble 
in the abrading environment represented by the coefficients 6, c, the box equa- 
tions suggest a natural generalization which tracks the mutual abrasion of two 
pebbles y, z, each of which represents the abrading environment for the other. 
(Alternatively, one may think of two pebble types abrading each other. We will 
return to the proper statistical interpretation of these equations in section |4]) 
Based on (fT2ll and (fT8|) we obtain the simultaneous ODEs: 

y = F(y,fe(z),c(z)) =F(y,z) (19) 
z = F(z,6(y),c(y)) =F(z,y) (20) 

describing the (approximate) evolution of box ratios under a purely collisional 
abrasion process. At first sight, the mutual abrasion of two pebbles may not 
seem to be a realistic model for the shape evolution of vast pebble collections 
on shingle beaches, nevertheless, we will show that equations p^ - (|20|) open up 
the possibility of describing the global interaction between abraded pebble and 
the abrading environment, including segregation. Another, important feature of 
the box equations is that they can be easily extended to model abrasion beyond 
collisional processes. First we look at two special collisional processes. 



2.1 The stationary case 

First we study the case where pebble z is infinitely harder then pebble y, i.e. 
we have z(i) = z(0) so ([T ^ - (|20p is reduced to the 3D system: 

y = F(y,z(0)), (21) 

which we call the steady state flow. Since here we have constant coefficients 
6, c, this case is directly comparable to the PDE ([3]). The infinitely hard pebble 
z could represent either the average value of a very hard pebble population 
in which the relatively soft pebble y is being abraded or, alternatively, the 
average value of a pebble population which remains invariant not because every 
pebble remains invariant but because it has a constant source and a constant 
sink. In the latter interpretation y is also a member of this population. We are 
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interested in the geometric properties of these flows. While the explicit formulae 
in (|T9l) - ([20|) admit a rather straightforward, rigourous mathematical study, here 
we merely state the main geometric facts as observations. 

As stated in section [51 the flow (PT|) can not have any genuine fixed points, 
and the only homothetic solution is the sphere S. Nevertheless, to explore the 
geometry of the flow we can still ask whether 'temporary' homothetic solutions, 
satisfying (I14[) may exist thus we seek the simultaneous solution of the quadratic 
equations 

ylFf + 2by,Ft' + cF^ = i = l,2. (22) 

Two out of four roots of (|22p are strictly positive, two strictly negative. The 
former define two surfaces Xi{yi, 1/2, b,c), {i — 1,2), where b,c are substituted 
from Above these surfaces (2/3 > A^) the flow moves away from S, towards 

0, below the surfaces (2/3 < Ai) the opposite happens, between the surfaces the 
derivatives 2/1, 2)2 have opposite signs. Points along the intersection lines defined 
by Ai — A2 satisfy ([Tl)) . trajectories hitting these lines slow down and turn 
back at cusps . Since the 1/2) projection of such trajectories loiters near the 
cusp, we call these points loitering points and the surfaces Ai, A2 we call loitering 
surfaces. The distance between the two loitering surfaces is rather small, so most 
trajectories turn at a sharp U-turn. Both loitering surfaces escape to infinity at 
O and both have coinciding, global minima at S: Ai(l, 1, &, c) = b + \Jb'^ + 3c. 
Based on this and by using ([5)). (|18l) we obtain a global, sufficient condition for 
the attractivity of the sphere: 

^ < zi + Z2 + 1 . (23) 

23 

In case of spherical objects with radius = R being abraded by spherical 
abraders with radius r {zi = 1, Z2 = 1,-23 = r) we arrive at the stability result 
of Bloore (|3]). In general, (I55|) is a global extension of this result (in the frame 
of the box equations); this condition does not depend on the box ratios 2/2 of 
the abraded object. In addition to the previous observations, one can also show 
that trajectories of pT|) transversally intersect each loitering surface Ai at most 
once. Combining these facts one can globally describe the steady state flow ((2T|) . 
The (2/1,2/2) projection of all trajectories of (f2T|) are either of 'Type F: O ^ S, 
or of 'Type 11': 5* — > 5, and they have at most one extremum in either variable, 

1. e. Type II trajectories are 'simple loops' originating and ending at the sphere 
(cf. Figure [5)/a). While the PDE ([3]) may have more complex structure, simple 
loops appear to be present. Figure [2]/b illustrates two trajectories originating 
at an ellipsoid with axis ratios j/i = 0.83, 2/2 = 0.91 evolved under ^ and (PT|) . 
respectively. While steady state flows correspond to coUisional pebble abrasion 
among well-mixed pebbles and this is an uncommon geological situation, there 
exist observations [33] which strongly suggest that non-monotonic evolution of 
box ratios may be present in some geological settings. A detailed description of 
this phenomenon is provided by Kuenen [19^ . 

We return to the experimental veriflcation of loops in section [3l We also 
remark that Type I. trajectories may correspond to the evolution of pebbles 
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Figure 2: Geometry of the steady state flows, a) Type I and Type II trajectories in 
the Zingg triangle. Turning points in j/1,1/2 are marked with L\,L2, respectively, b) 
Simulation of Type I trajectory: ellipsoid with axis ratios y\ = 0.278, 3/2 = 0.329, j/3 = 
1.00 evolved under the PDE (|3]) and under the steady state box flows ^21\ in the 
presence of constant spherical abraders with radius r = 0.5 {zi — Z2 — 1, — 0.5, b = 
0.5, c = 0.25) c) Simulation of Type II trajectory: ellipsoid with axis ratios yi = 
0.83,2/2 = 0.91, j/3 — 1.00 evolved under the PDE (|3]) and under the steady state box 
flows (|2ip in the presence of constant spherical abraders with radius r = 0.1 



carried by rivers where the abrading environment is the riverbed, represented 
by a very hard and very large abrader. 

2.2 The self-dual case 

Now we study the case where two identical pebbles y, z mutually abrade each 
other, i.e. we have y{t) = z{t) so (fT9 |) -([20 |) is again reduced to a 3D system: 

y-F(y,y), (24) 

which we call the self-dual flow with no free parameters and an autonomous 
direction field in the (2/1,2/2) plane: 

where fy are defined in pS]) by a z = y substitution and the superscript C 
stands for 'collisional'. Since (|24p has time-dependent coefficients, it is difficult 
to compare it directly to Nevertheless, it represents the abrasion of pebbles 
in an environment similar to itself, i.e. it is the simplest model of the abrasion 
of segregated pebble populations. The geometry of self-dual flows is simple 
(visually rather similar to the Gaussian flows illustrated in Figure [ij , in fact all 
trajectories are of Type I and converge globally to the sphere and in we 
have 

F,^y,>0 i^l,2 (26) 

on all trajectories, at all times. This is not very surprising, not only because 
it follows from (|25p but also because it intuitively matches The latter 
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condition was derived for constant z, however, we can see that it is fulfilled in 
the self-dual flow on every trajectory, at any time. Since self-dual flows are the 
simplest model including segregation, they are closer to the physical process 
of pebble abrasion on shingle beaches than stationary flows. While self-dual 
flows have an autonomous 2D component as in (ITS]) we showed that such a 
system might produce homothetic solutions, nevertheless in the collisional box 
equations there is no such indication. 



3 Friction 

3.1 Basic assumptions and a semi-local PDE model 

From the previous sections we learnt that collisional abrasion, as modeled by 
eq. (O is dominated by the Gaussian and Mean Curvature flows (|10p thus it 
converges ultimately to the sphere. It is not entirely surprising that collisional 
abrasion is determined by local quantities since in this model we assume that 
abrading particles arrive from uniformly distributed directions. This model is 
based on the assumption that in the abrasion process the energy level is high 
enough to support the free flight of mutually abrading particles. 

Frictional abrasion is rather different: here the kinetic energy can be rather 
low and the orientation (and thus the global shape) of the particle play a key 
role so we look for a semi-local PDE model depending both on local and on 
global quantities. Adopting Aristotle's approach, we regard the radial distance 
R (measured from the ccntroid) as our primary variable but we also assume fric- 
tional abrasion to depend on the global minimum and maximum of R described 
by the variables Rmim Rmax, respectively: 

g'^ ' ' ^ ~f{R,R„,^n,RmaxMu,V,t) , />0. (27) 

Thus our model may be regarded as a generalization of that of Aristotle ([T]). 
Friction influences abrasion either by sliding or by rolling, and we proceed by 
making specific assumptions about these components, starting with sliding. In 
shear band experiments with granular materials it has been observed ( 31 ) that 
the long axis of elongated particles gets aligned with the direction of shear. This 
result indicates that in case of sliding, abrasion will be concentrated on the fiat 
sides of the abraded pebble. This implies that for sliding friction we can assume 

In case of rolling abrasion the situation is radically different. Rolling is essen- 
tially a 2-dimensional motion, the plane of which is determined by the kinetic 
energy. If this energy is so high that the pebble can even roll over the point 
furthest from the centroid then we are at the transition to collisional abrasion. 
There is also a lower energy threshold below which rolling is not possible: this 
corresponds to a trajectory with minimal maximum distance. Rolling occurring 
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between these two critical energy levels results in a non-monotonic abrasion law, 
i.e. we assume that abrasion speed is small for radii to the minimum, zero at 
the maximum and larger in-between, i.e we assume the existence of a radius R^,, 
Rmin l£ R* l£ Rmax for which wc havc 

< /(i?™n) < I{R.) > f{Rma.) - . (29) 

The non-monotonicity of the abrasion law could be explained by observing that 
surface points corresponding to very large radii are hardly affected, because they 
are not part of the typical trajectory. Points corresponding to very small radii 
are part of typical trajectories, however, they will abrade less than points with 
larger radii along the same trajectories. In other words: rolling is essentially 
a planar problem, the plane of rolling prefers smaller radii because of energy 
minimization, however, from among radii in the plane of rolling, larger radii will 
be abraded faster (in accordance with Firey's ([IH]) argument). 

The above ideas can be translated into what we shall call a semi-local PDE 
model, i.e. we propose a definite form for for / in (j27l) . Here we introduce 
separate terms for sliding and rolling with independent coefficients Vsji^n re- 
spectively and use the dimensionless ratios ri = -R/i?mm,^2 = R/Rmax- 

f{R, Rmin, Rmax) Vsfs{ri,r2) + l^rfr{ri,r2) = I^s?'27'r" + '^rr2il - ^2 ) • (30) 

For sufficiently high values of n, this model appears to capture most essential 
physical features of the processes we are aiming to describe. The PDE model is 
clearly not unique but provides a simple basis for a qualitative analysis. 

3.2 Unified model for collisional and frictional box flows 

In the box approximations, assumptions ()28p and (|29p translate into 



ui < U2 < W3 < 0, U2 l£ ui < ii^ < 0, (31) 

respectively. In the n — ^ oo limit the PDE defined by (|27p . (l30l) yields for the 
box flows 

Ul = -VsUl - VrUl, U2 = -Vry2, ^3 = , (32) 

which is equivalent to 

y^Ff{y,i^„iyr) = -iiysF^ + i^rF'^), (33) 

where 

F^ = -[yi,0,0]^ , F^ = -[2/i,y2,0]^. (34) 

These equations satisfy the conditions (|3T|) . Similarly to the self-dual flows 
in the collisional model, we can observe that while the full flow is 3D, both 
component flows (Sliding, Rolling) can be reduced to simple planar flows: 
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and in p3p we have 

F/ = y,<0 z = l,2. (36) 

We can now simply add collisional and frictional flows (IT^ - (P0)) . to obtain 
the unified equations (denoted by superscript u) for the two-body problem: 

y = F"(y,z,i.,,M,) =F(y,z)+F-^(y,i.„i.,) (37) 
i = F"(z,u,i/„iy,) =F(z,y)+F^(z,zy„i/,). (38) 

3.3 The steady-state flow in the presence of friction: ex- 
perimental evidence of simple loops 

We study the steady-state flows in the presence of friction, defined by the 3D 
system: 

y = F"(y,z(0),i/„i^,) =F(y,z(0))+F^(y,i.„z/,). (39) 

The basic geometry is similar to the collisional case, described in subsection 
12.11 we can only observe Type I {O ^ S) and Type II {S — S) trajectories. 
To verify the latter we conducted a simple table-top experiment in the spirit of 
Kuenen jl7j. In each experiment we rolled 5 rounded chalk pieces of approxi- 
mately 9.8mm diameter surrounded by approximately 400 hard, plastic balls of 
diameter 1.8mm in a horizontal glass cylinder at 5 rpm for approximately 24hrs. 
The box dimensions ui,U2,U3 were measured every hour, {u^ was the largest 
diameter, U2 the largest diameter orthogonal to M3 and ui was the diameter 
orthogonal both to and U2.) The data for the 5 chalk pieces was averaged. 
For the computations we used the unified frictional equations p9p and the co- 
efficients b,c were computed from ^TE\\ . In case of abrading plastic balls with 
radius r — 0.9 we had b = r — 0.9, c — r'^ = 0.81, respectively. Initial conditions 
were set identical to the experiments and the friction coefficients vs,vr were 
selected to achieve a good match. 

As we can observe in figure |3l both experiments and numerical simulations 
show fair quantitative agreement, both displaying the simple loop predicted in 
subsection 12.11 



3.4 The self-dual flow in the presence of friction and the 
existence of stable, nontrivial attractors 

Analogously to the self-dual collisional fiows (|24l) we define the self-dual fric- 
tional flows as the mutual abrasion of two identical pebbles in the presence of 
friction: 

y = F"(y,y,i.,,i.,) =F(y,y)-|-F/(y,i.„z.O- (40) 
which, similarly to (j25p . also have an autonomous 2D direction field, however, 
with two free parameters Vs,i^r- We look for homothetic solutions ([T3|): as fixed 
points of (PUI) . i.e. we solve ~ F2 = which is a 2x2 linear system for Vs^Vrj 
yielding: 

_ FiFi'-F2Ff F2Fl'-F,F^ 



^1^2^ — F2Ff ' F1F2' — F2F^ 



(41) 
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EXPERIMENTS 

BOX SIMULATIONS 






Inltiat conditions: 
y_l (0)=O.9B3 
V_2(D)=a.988 




Friction coelTiciemts: 

v_s=0.21 

v_i=0-63 




Time step: 
t«=D.00O5 




Time: 

11 = 0.432 


yi 



0.6 0.7 0.8 0,9 1 



Figure 3: The stationary flow in the presence of friction. Observe the marked 
'simple loop' both in the numerical and the experimental data 



where Fi are from ([8|), Ff,Ff^ are from ([34|) . By substituting (yi = yio,2/2 = 
j/2o) into (HI])) we obtain the friction coefficients in the presence of which (|40p 
will have a homothetic solution at (yi — yw,y2 — y2o)- 




Figure 4: 'Bracketing' conditions for the existence of homothetic solutions, a) 
Equilibrium in the collisional flows, b) Equilibrium in self dual flow in the 
presence of friction. 
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Friction will produce physically observable, stable attractors if both critical 
parameters are positive. The condition for this is analogous to the necessary 
condition (jlSp defined for purely collisional flows. Here we regard the vector sum 
of the self-dual collisional flow with the frictional flows. To obtain equilibrium, 
the former must be 'bracketed' by the latter; this can be expressed based on 
1^ and (133 as (cf Figure [l)/b): 

-yf)(yf -2/f)<0. (42) 

Condition (|^^ is satisfied by our equations (|33I) - ([M)) for friction, in fact, the 
stricter condition 

= yf <yf <^<yf (43) 

yi 

is also met, and therefore the proposed PDE ([50]) is capable of predicting the 
existence of nontrivial attractors in the (j/i, 2/2) plane. Since both Fi, F2 appear 
to be monotonically decreasing in yi,y2, the fixed points defined by (j4ip are 
stable attractors, illustrated in Figure [5]/a. 




Figure 5: The self-dual flow in the presence of friction with nontrivial critical point 
at (2/1,2/2) ~ (0.5,0.8). a. Friction law according to (|30|l : single, stable nontrivial 
attractor, parameter values {vs ~ 0.622, Ur ~ 1.909). b. Friction law according to 
l|47[l: nontrivial saddle, both the sphere at (2/1,2/2) = (1, 1) and infinitely flat shapes 
at 2/1 = are attractive. Parameter values {us — 1.244, i/,. — 1.909). 



4 Stochastic process: multi-body simulations 

The unified collisional- frictional equations ([57|) -([55 1) describe the mutual abra- 
sion of two individual pebbles, however, they also define a Markov process where 
y, z are random vectors with identical distributions since they represent two 
random samples of the same pebble population. The evolution of this Markov 
process (and thus the time evolution of of the pebble size and ratio distribu- 
tions) is of prime interest since it determines the physical relevance of the stable 
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attractors identified in subsection 13.41 While the analytical investigation of the 
Markov process is beyond the scope of this paper, direct simulations are rela- 
tively straightforward. We consider N pebbles out of which we randomly draw 
two with coordinates y°,z° and run equations (|57)) - (|55)) for a short time pe- 
riod At on these initial conditions to obtain the updated vectors y^,z^ . In the 
simplest linear approximation we have the recursive formula 



Such an iterative step can be regarded as the cumulative, averaged effect of sev- 
eral collisions between the two selected pebbles. Apparently, the iV = 2, Ai -> 
case is identical to ([37l) -(|38 |) . Previously we investigated the behaviour of the 
deterministic flows in the special cases of steady state and self-dual flows. Multi- 
body simulations allow the numerical study of the statistical stability of the 
flows, i.e. one can assess the stability of the above-mentioned special cases. 

4.1 Stochastic, multi-body simulation of collisional flows 

The characteristic "U"' turn identified in the continuous flow for Type II tra- 
jectories appears to be very robust. In multi-body simulations one can observe 
the fuzzy zig-zag geometry of the trajectory near the turning point indicating 
that larger amount of time is spent near the 'loitering' points. This behaviour 
is illustrated in Figure ^a. for N=4. The self-dual collisional flow proves 
to be remarkably stable and in multi-body simulations the trajectories of the 
continuous flow are often followed rather closely. However, there are exceptions. 
Figure |6j^b illustrates an N=9 body simulation with random initial conditions, 
where 8 trajectories converge to sphere (similar to the continuous flow), how- 
ever, one trajectory turns back. This deviation from the continuous case can be 
expected in the vicinity of the sphere since here the derivatives yi , 1/2 converge 
to zero and so small perturbations can change their signs. 

4.2 Stochastic, multi-body simulations of the unified col- 
lisional and frictional model 

In subsection 13.41 we showed that in the presence of friction stable, nontrivial 
attractors appear in the self-dual flows (|40l) . While these fix points collectively 
attract the abrading objects, the latter mutually repel each other as we illus- 
trate in Figure [71 where two almost-identical pebbles initially follow the same 
trajectory, however at some critical time they split. Here we used (P^ - (|^5|) with 
N — 2 and almost identical initial conditions to simulate (1401) . 

This shows that self dual fiows (describing the mutual abrasion of a ho- 
mogeneous pebble population consisting of identical pebbles) are structurally 
unstable in the sense that the homogeneity of initial data decays (pebble ratios 
may mutually repel each other in the abrasion process). Based on this we expect 
that in multi-body random simulations according to p4)) - (|45|) the attractors do 




AtF"(y',z') 



(44) 
(45) 
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Figure 6: Randomized simulation of collisional flows, a) 4-body simulations of 
the trajectory's 'U' turn in the steady state flow, b) 9-body simulation of the 
self-dual flow with random initial conditions. 8 trajectories approach the sphere 
(as in the deterministic case) , one trajectory turns back (see right hand enlarged 
version) . 



not appear. The geological motivation for the study of self-dual flows was the 
observation that pebbles are often segregated by size by the wave transport. If 
this process is happening on the same (or faster) timescale as abrasion, then 
the self-dual flows and their attractors may be stabilized by sustaining the ho- 
mogeneity of the initial data, i.e. by selecting only nearly similar pebbles for 
the collision process. This could be implemented by introducing a correlation r 
between the random vectors y, z; our observation of the instability corresponds 
to the uncorrelated r = case. The fully correlated r = 1 case is essentially 
identical to the deterministic evolution of ()40p . Another option to implement 
segregation (possibly closer to the geological process) is to omit pebbles from 
the simulation the maximal size y^^i of which differs from all other pebble sizes 
y3,j by a prescribed factor. For each pebble the coefficient 



1 I y3.i 
log — 



(46) 



was computed and pebbles with (Ji > a were omitted. Increasing segregation 
effect corresponds to decreasing values of cr; this is illustrated in Figure [51 Ob- 
serve that for weak segregation a — 2 the stable fixed point of the deterministic 
flow vanishes while for stronger segregation it persists. 
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Figure 7: Splitting of trajectory in the stochastic simulation (|44 |) - (|45|) of unified 
collisional and frictional fiows (HUl) , created by the mutual abrasion of two almost 
identical particles. Critical point of the deterministic fiow (cf. Figure [5] ) at 
(2/1,2/2) — (0.5,0.8) at parameter values {vg — 0.622, — 1.909)). 




Figure 8: The self-dual fiow in the presence of friction and segregation. Critical 
point of the deterministic flow (cf. Figure[5]) at (2/1,2/2) = (0.5,0.8) at param- 
eter values {vs = 0.622, Vr = 1.909) For weak segregation (tr = 2) the random 
trajectories do not converge to the fixed point, however, with increasing segre- 
gation the fixed point of the deterministic fiow becomes attractive also in the 
stochastic process 



5 Conclusions 

Our goal was to create a framework in which the interaction of mutual particle 
abrasion with global transport and segregation can be studied. 
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5.1 Mutual abrasion versus individual abrasion 

As a first step we radically simplified Bloore's classical PDE, describing the 
general collisional abrasion process of an individual pebble. Our system of 
ODEs, called the box equations, attempts to track the evolution of maximal 
pebble size and main pebble proportions. While the box equations are merely a 
heuristic simplification of the PDE describing individual abrasion, they admit 
the global, analytic study of mutual abrasion. To achieve this, we incorporated 
the results of Schneider and Weil [32], Varkonyi and Domokos providing 
the general interpretation of the coefficients in Bloore's PDE. We arrived at the 
system of 6 coupled ODE's , defining the mutual abrasion of two pebbles. We 
investigated two special cases of this rich process when the 6D system reduces 
to a 3D system: the steady-state flows where one pebble is regarded as constant, 
and the self-dual flows where both pebbles and abraders are treated identically. 

5.2 Special cases: steady state and self-dual flows 

Steady state flows may serve as a skeleton model of abrasion of well mixed pebble 
populations where the abrading environment has stationary size distributions. 
In the box equations we extended Bloore's local stability result for nearly- 
spherical objects abraded by spherical objects to obtain the sphere's global 
attractivity criterion. Equation ()23|) provides a sufficient, global criterion for 
the convergence to the sphere, depending only on the size ratio of the abraded 
and abrading object and the box proportions of the latter. We found that the 
steady state box flows admit only two types of trajectories, one of which appear 
as simple loops in the plane of the two box proportions, originating and ending 
at the sphere. We confirmed the existence of this 'boomerang effect' by direct 
simulation of the PDE, in simple table-top experiments and also in stochastic 
simulations. While loops may cause some statistical accumulation of pebble 
ratios in the vicinity of the turning point (which we called 'loitering points'), 
they certainly do not account for the very clear and striking appearance of 
dominant pebble ratios. 

We found that the self-dual flows, modeling the abrasion of pebbles in an 
environment identical to the abrading pebble, uniformly and globally converge 
to the sphere. These flows represent the box equation model for pebble abrasion 
in segregated environments. 

5.3 Friction 

Since none of the collisional equations accounted for the existence of homothetic 
solutions, we added frictional terms, corresponding to sliding and rolling, with 
coefficients lys^i^r, respectively. We postulated a semi- local PDE (pO)) based on 
simple physical observations and derived box equations describing the unified, 
self-dual collisional and frictional process. We found that even arbitrarily small 
amount of friction can stabilize nontrivial shapes as homothetic solutions in the 
vicinity of the sphere. We also observed that sliding friction stabilizes discoid 
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shapes whereas rolling friction stabilizes elongated rod-like shapes. At larger 
values of the coefficients Vs^i^r, shapes further away from the sphere can be 
stabilized. Physically, the increase of the coefficients might be due either to 
the increased importance of frictional efficiency or, more importantly, to the 
decrease of the total energy of the abrasion process. In the box equations we 
determined the coefficients Vg , corresponding to fixed points with given lo- 
cations. The latter data could be extracted from statistical measurements on 
shingle beaches and our equations could be used to identify the associated coeffi- 
cients. Our proposed semi-local PDE model is certainly not unique, alternative 
PDE's , also meeting the same criteria, can predict multiple attractors for the 
same value of the coefficients; therefore several dominant box ratios could co- 
exist in a pebble environment governed by mutual abrasion by collisions and 
friction. As an example we mention the following model as an alternative to 

/(i?,i?™„,i?„a.) = ^^.i^rr"-)-z.,ri^^l^(l-r^), (47) 

r2 ri 

which, similarly to (1501) . meets the criteria (1^51) and (P^ : the corresponding box 
flow (defined in the n — > oo limit) is illustrated in Figure [n]/b with a saddle- type 
critical point. 



5.4 Statistical approach: segregation by transport 

While the existence of attractors in the presence of friction is promising, the 
model calls for statistical verification. Our box equations in the self-dual case 
describe the mutual abrasion of two identical pebbles, i.e. the abasion process 
where the abraded pebble and abraders are treated completely symmetrically. 
To analyze the statistical stability of the results, we described the Markov pro- 
cess based on the box equations. This process defines the evolution of pebble 
size and ratio distributions under the combined action of collisional and fric- 
tional abrasion and segregation. We ran direct simulations starting with almost 
identical initial conditions. In this process, the colliding pair of pebbles is drawn 
randomly from the same distribution. First we modeled abrasion without trans- 
port and segregation. To this end, the draw was uncorrelated and we found that 
the attractors disappeared: this matched the geological observations of Carr [6], 
Landon |20j and Bluck [5, claiming that global transport played an important 
role in the emergence of equilibrium shapes. Next we introduced transport and 
segregation by size, so the maximal size of the colliding pairs was highly corre- 
lated and we found that stable attractors were sustained in the random process. 
Our model thus predicts that dominant pebble ratios as stationary equilibrium 
shapes may robustly appear in abrasion processes where, in addition to colli- 
sions, friction also plays an important role, and an ongoing global transport 
process segregates pebbles based on size. Segregation by shape (also included 
in our model) would enhance the stability of these dominant ratios even fur- 
ther. These predictions match well geological observations, also, they confirm 
Aristotle's intuition that the basic mechanisms of pebble abrasion can only be 
understood by including semi-global and global effects. 
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